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ABSTRACT 

We analyse archival Rossi X-Ray Timing Explorer (RXTE) proportional 
counter array (PCA) data of thermonuclear X-ray bursts from the 2002 outburst 
of the accreting millisecond pulsar SAX J1808. 4-3658. We present evidence of 
a complex frequency modulation of oscillations during burst rise, and correla- 
tions among the time evolution of the oscillation frequency, amplitude, and the 
inferred burning region area. We discuss these findings in the context of a model, 
based on thermonuclear flame spreading on the neutron star surface, that can 
qualitatively explain these features. From our model, we infer that for the 2002 
Oct. 15 thermonuclear burst, the ignition likely occured in the mid-latitudes, the 
burning region took ~ 0.2 s to nearly encircle the equatorial region of the neutron 
star, and after that the lower amplitude oscillation originated from the remaining 
asymmetry of the burning front in the same hemisphere where the burst ignited. 

We emphasize that studies of the evolution of burst oscillation properties dur- 
ing burst rise can provide a powerful tool to understand thermonuclear flame 
spreading on neutron star surfaces under extreme physical conditions. 

Subject headings: equation of state — methods: data analysis — stars: neutron 
— X-rays: binaries — X-rays: bursts — X-rays: individual (SAX J1808. 4-3658) 


1. Introduction 

X-ray bursts are produced by thermonuclear burning of matter accumulated on the 
surfaces of accreting neutron stars (Grindlay et al. 1976; Belian, Conner, & Evans 1976; 
Woosley, &; Taam 1976, Joss 1977; Lamb, & Lamb 1978). During many bursts, millisecond 
period brightness oscillations are generated by the combination of rapid stellar rotation 
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and an asymmetric brightness pattern on the neutron star surface (Strohmayer, k Bildsten 
2003). The period of these oscillations is very close to the stellar spin period (Chakrabarty 
et al. 2003; Strohmayer et al. 2003). Moreover, as this timing feature originates from the 
surface of the neutron star, its detailed modeling may be useful to constrain stellar structure 
parameters, and hence the equation of state models of the dense matter in the neutron star 
core (Bhattacharyya et al. 2005; Miller, k Lamb 1998; Nath, Strohmayer, k Swank 2002; 
Muno, Ozel, k Chakrabarty 2002). 

Modeling of burst oscillations can also be useful to understand neutron star atmospheres, 
surface fluid motions, and for mapping the magnetic field structure on the stellar surface. 
For example, the evolution of frequency of these oscillations during the burst rise phase, may 
provide information on the spreading of thermonuclear flames under the extreme physical 
conditions that exist on neutron stars (Bhattacharyya k Strohmayer 2005b). This is because 
bursts almost certainly ignite at a particular point on the stellar surface (as simultaneous 
ignition over the whole surface would require very fine tuning), and then spread to burn 
all the surface fuel (Spitkovsky, Levin, k Ushomirsky 2002; Bhattacharyya k Strohmayer 
2005a). This slow (compared to the rotational speed) movement and spreading of the burning 
region may give rise to complex frequency evolution of the observed burst oscillations. This 
spreading would also cause the observed burst intensity to increase, and the oscillation 
amplitude to decrease. Moreover, the increase in emission area can be estimated by spectral 
analysis (Strohmayer, Zhang k Swank 1997; Bhattacharyya k Strohmayer 2005a; 2005c). 
Therefore, simultaneous modeling of the evolution of burst intensity, oscillation frequency 
and amplitude, and spectral properties can, in principle, be a powerful tool to understand the 
propagation of burning fronts on neutron star surfaces under conditions of extreme radiative 
pressure, magnetic field and gravity. 

Frequency evolution during burst rise oscillations has so far been observed from two low 
mass X-ray binary (LMXB) sytems: SAX J1808. 4-3658 and 4U 1636-536 (Chakrabarty et 
al. 2003; Bhattacharyya k Strohmayer 2005b). The Rossi X-ray Timing Explorer (RXTE) 
observed the 401 Hz X-ray pulsar SAX J1808. 4-3658 in October and November of 2002, 
when it was in outburst. Four type I X-ray bursts were detected during these observations 
(Chakrabarty et al. 2003), three of which showed strong millisecond period brightness os- 
cillations during burst rise. A previous study found that as the burst intensity rises, the 
oscillation frequency also increases by ~ 5 Hz and may overshoot the stellar spin frequency 
(Chakrabarty et al. 2003). Here, we analyse these archival data in order to model the time 
evolution of different burst properties and search for correlations among them. 

In our study we find evidence of complex variations in the oscillation frequency during 
the rising phase of bursts. This is the first report of such variations from any source. The 
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frequency modulation is correlated with the evolution of oscillation amplitude and burning 
region area. In § 3, we propose a qualitative model, based on thermonuclear flame propaga- 
tion, that can explain the evolution of different burst properties, as well as the correlations 
among them. 


2. Data Analysis and Results 

We analyse the archival RXTE proportional counter array (PCA) data of the 2002 
outburst from SAX J 1808. 4-3658. Three thermonuclear bursts with significant millisecond 
oscillations during the rising phase are found in the Obslds: 70080-01-01-000 (Oct 15), 70080- 
01-02-000 (Oct 18), and 70080-01-02-04 (Oct 19). First we explore the frequency evolution 
of burst rise oscillations during these bursts using three procedures. We calculate dynamic 
power spectra (Strohmayer & Markwardt 1999) with time sampling short enough to resolve 
the burst rise interval, but large enough to accumulate sufficient signal power above the 
noise level. The dynamic spectra (panel a of Fig. 1, and Fig. 2) provide indications of 
the frequency evolution behavior. To confirm the indications in the dynamic spectra, we 
carry out a phase timing analysis (Muno et al. 2000). We divide the burst rise time interval 
into several bins of a fixed chosen length, and then assuming a frequency evolution model, 
we calculate the average phase (-0*,) in each bin ( k ). The corresponding x 2 is calculated 
using the formula x 2 = J2b=x(i>k ~ ipk) 2 /<4 (Strohmayer & Markwardt 2002), where M 
is the number of bins. For this study we used extensive burst rise simulations to evaluate 
the uncertainty, a^ k , as a function of the Z 2 power in each time bin. We find the best fit 
parameter values for various frequency evolution models by minimizing y 2 , and we calculate 
the uncertainty in each parameter by finding the change which produces the appropriate 
increase in x 2 (Strohmayer & Markwardt 2002; Press et al. 1992). Finally, we calculate 
the total Z 2 power (Strohmayer & Markwardt 2002) for the entire rise interval (ie. without 
binning) using the best fit frequency evolution model parameters, and ensure that this power 
is close to the maximum power obtained from any parameter values. 

We first fit the oscillations during the rising phase of the Oct. 15 burst with a constant 
frequency model. This gives a best fit frequency of 400.91 Hz and x 2 /dof = 30.04/8. We 
evaluate the significance of this x 2 value using simulations, and find a probability of 0.005 
to obtain such a value by chance. Since the constant frequency model does not describe the 
data well, we next consider more complex models. These are; (1) linear frequency increase, 
(2) second order polynomial, and (3) linear increase and subsequent linear decrease models 
(see Table 1 for a description of the models). These models give x 2 /dof values of 17.73/7, 
17.35/6, and 11.68/5 respectively. Although these models give better fits, they still have 
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reduced x 2 values > 2, and do not describe the data very well. A slightly more complex 
model which fits the data well has a linear frequency increase followed by a second order 
polynomial, and gives a x 2 /dof value of 3.36/4. Of the models tested we consider this the 
best description of the Oct. 15 burst, and the best fit parameter values are given in Table 
1. Panel a of Fig. 1 shows that this model is consistent with the dynamic power contours, 
and it indicates a frequency increase (by a few Hertz) for the first ~ 0.2 s from burst onset, 
then a frequency decrease (by ~ 1 Hz), and a subsequent increase. From Table 1, we note 
that the model parameters Vq & *d are practically unconstrained from the lower and upper 
sides respectively. This implies that the initial frequency increase can be well fit by a very 
steep model. However, the other sides of these parameters are reasonably well constrained, 
which shows that a constant frequency model is insufficient to model this portion. 

Next we fit the rise oscillations from the Oct. 18 & Oct. 19 bursts with the same 
frequency evolution models. The constant frequency model for the Oct. 18 burst gives a 
X 2 /dof = 154.51/8, and hence can be strongly rejected. This confirms the conclusions of 
Chakrabarty et al. (2003) who first noted the large frequency increase present during this 
burst. A linear frequency increase model gives a x 2 /dof = 19.35/7 that, though better, is still 
uncomfortably large to be acceptable. The next more complex model (a constant frequency, 
& subsequent linear increase) gives a x 2 /dof = 8.44/6 (see Table 1, and the upper panel of 
Fig. 2), and is statistically acceptable. However, we note that the dynamic power contours 
(upper panel, Fig. 2) are suggestive of an initial frequency increase, decrease, and increase 
behavior (shown by the dotted curve) qualitatively similar to that for the Oct. 15 burst. 
Although higher signal to noise ratios per time bin would likely be required to confirm this 
behavior, this (dotted curve) model does give a higher total Z 2 power than that given by 
the best fit x 2 model (solid curve) in Table 1. Finally, the constant frequency model for the 
Oct. 19 burst gives a x 2 /dof = 22.03/8, which is also unacceptably high. A linear frequency 
increase model (see Table 1, and the lower panel of Fig. 2), gives a x 2 /dof = 6.46/7, and is 
acceptable for this burst. 

Panel b of Fig. 1 shows the rms amplitude variation with time during the rise of the Oct. 
15 burst. From the beginning of the burst, the amplitude decreases for ~ 0.2 s, and then 
assumes a near constant value, with some fluctuations. This behavior is qualitatively similar 
to that seen in bursts from the LMXB systems 4U 1728-34 and 4U 1636-536 (Strohmayer, 
Zhang & Swank 1997; Bhattacharyya & Strohmayer 2005b). We also perform time resolved 
spectral fitting (using a blackbody model) during the rise of the Oct. 15 burst. The inferred 
source radius (which provides some relative indication of the size of the burning region on the 
stellar surface) shows evidence for a modest increase for the first ~ 0.2 s, and then remains 
almost constant (panel c, Fig. 1). Panel d gives the corresponding source temperature 
variation. 
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3. Discussion 

Taken together our results for the three bursts from SAX J1808. 4-3658 suggest that 
the oscillation frequency can evolve in a complex manner during burst rise. In all cases a 
constant frequency model is a poor description of the oscillations. In two bursts (Oct. 15 
and 18) the data can best be described by a complex modulation in frequency whereby the 
frequency initially increases, then decreases before increasing again. We now discuss from 
a theoretical perspective how spreading of thermonuclear flames can plausibly account for 
such complex modulations in the oscillation frequency and amplitude. In our discussions we 
use the observational results mostly from the Oct. 15 burst as characteristic. The salient 
features of this burst are as follows: (1) from the start of the burst, the oscillation frequency 
and burning region area increase and the oscillation amplitude decreases for ~ 0.2 s; (2) after 
~ 0.2 s, the frequency first decreases and then increases, and both amplitude and burning 
region area reach a nearly constant value (with some fluctuations) . 

The burst begins when the fuel (i.e., accumulated matter) ignites at a particular point, 
and then the flame propagates over the surface (Spitkovsky et al. 2002; Bhattacharyya 
& Strohmayer 2005a; 2005c). Before spreading has engulfed the entire star, temperature 
variations due to surface waves may not be able to explain the brightness oscillation or its 
frequency evolution (as has been proposed for the burst tail oscillations; Heyl 2005; Lee & 
Strohmayer 2005), as the rapid spreading and temperature increase may wash out this effect. 
This explanation of frequency evolution during burst rise was also shown to be unfavored 
by Bhattacharyya & Strohmayer (2005b). However, thermonuclear flame propagation can 
give rise to the observed frequency evolution, as different accelerations of the eastbound and 
the westbound burning fronts can cause acceleration (either eastwards or westwards) of the 
center of the burning region relative to the star. This is our primary explanation for the 
observed frequency modulation. Here, we do not consider oscillations which could be caused 
by local vortices or magnetic fields, as even if they form, they may not necessarily generate 
a hot spot. 

We first review some relevant results from previous work (Spitkovsky et al. 2002): (1) 
the greater scale height of the burning region (than the cold fuel) gives rise to a shearing 
speed (as the horizontal pressure gradient in the burning front increases with height). As a 
result, the cold fuel is drawn into the burning front and ignited. This enables the burning 
front to propagate. (2) The shearing speed is greater nearer the equator than the pole (due 
to the latitude dependence of the Coriolis parameter; Spitkovsky et al. 2002). Thus, the 
burning front propagates with the shearing speed ($) (assuming the mixing time scale is 
very small; Spitkovsky et al. 2002; Fujimoto 1988; Cumming & Bildsten 2000). Now, if 
the fuel ignites at a mid-latitude (say, in the northern hemisphere), it will propagate faster 
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towards the equator (than the pole), and the east- west width of the burning region will 
increase much faster near the equator. Therefore, after a certain time, the burning region 
will encircle the equator and propagate more or less symmetrically towards the south pole. 
The northern burning front propagates towards the north pole, keeping the asymmetry (due 
to the variation of east- west width with latitude), which vanishes near the pole. At the 
beginning, the burning region is relatively small, and hence the oscillation amplitude can be 
large. As the burning region grows, the oscillation amplitude naturally diminishes (panel b 
of Fig. 1). This effect was also reported for bursts from other sources (Strohmayer, Zhang, 
& Swank 1997). After the burning region encircles the equator with a considerable north- 
south width, the observed burning area does not increase much (hence the near constant 
radius after the initial increase; panel c of Fig. 1). From this time, the oscillation is due 
to the asymmetry of the northern burning front (with the persistent background due to the 
azimuthally symmetric portion of the burning region), and hence the amplitude attains a 
near constant value till the asymmetry vanishes. Therefore, according to our explanation, 
for the Oct. 15 burst, the burning region takes ~ 0.2 s to nearly encircle the stellar equator. 

We consider two additional effects to explain the complex frequency evolution: (3) due 
to the increased scale height, the top portion of the burning region slips westwards (if the star 
rotates eastwards) to conserve angular momentum. If the shearing speed due to this is v, the 
burning front propagates westwards with a speed v+tf and eastwards with a speed d. (4) The 
angular momentum of a stellar surface element increases from pole to equator. Therefore, to 
conserve angular momentum, if the burning front propagates towards the equator, it drifts 
westward, and if it propagates towards the pole, it drifts eastward. The corresponding speed 
( V ) depends on latitude. For a given value of d, V is larger near the pole, as the fractional 
change of angular momentum is larger there for a given amount of latitude change. However, 
V is a monotonically increasing function of d, and d increases from pole to equator. So, 
two competing effects determine the latitude and direction dependence (northbound and 
southbound) of V. Another effect may also be relevant for the direction dependence of V. 
At a given latitude in the northern hemisphere, if the front moves southwards, it is pulled 
downwards by the increase in the effective gravity (as it conserves its angular momentum), 
perhaps making the friction and the mixing within the front more efficient. Hence the front 
speed both southwards and westwards may increase, and the value of V is greater. For the 
opposite reason, at the same latitude, a northbound front has a smaller V . Considering 
these three effects on V, it seems likely that in the northern hemisphere, V is larger for a 
southbound front at low-latitudes than for a northbound front at mid-latitudes. 

Our explanation for the observed frequency modulation uses all four effects discussed 
above. If ignition starts at a mid-latitude in the northern hemisphere and initially propagates 
mostly southwards, its eastbound front has a speed while the westbound front has speed 
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V + v + (the star rotates eastwards). Therefore, as $ is likely to be smaller than V + v 
in the mid-latitudes, the observed frequency starts from a value ~ a few Hz lower than the 
stellar spin frequency. As the front approaches the equator, increases quickly, and as it 
crosses the equator, the term V switches from westbound to eastbound. These two effects 
quickly bring the observed frequency back towards the stellar spin frequency. Indeed, if V 
is large enough then it may cause the observed frequency to overshoot the spin frequency. 
This goes on till the burning region more or less encircles the star around the equator. The 
fact that the peak frequency occurs ~ 0.2 s (panel a of Fig. 1) after the beginning of the 
burst supports this qualitative scenario. After the equatorial belt has ignited, the oscillation 
results from the remaining asymmetry of the northern front. The westward speed of the 
west front of this asymmetry is v + while the eastward speed of its east front is V + d. 
The competition between these two speeds determines the frequency evolution. If initially, 
during the time interval when this asymmetry takes over, v > V (in the mid-latitudes), then 
the observed frequency decreases, and then, as V eventually dominates (as the front proceeds 
northwards), the observed frequency starts to increase again, and may, in fact, overshoot 
the stellar spin frequency. This can explain the frequency evolution of the Oct. 15 burst 
(panel a, Fig. 1). However, if V > v immediately after the equatorial belt is ignited, the 
oscillation frequency increases monotonically (lower panel of Fig. 2; also Fig. 1 & Fig. 2 
of Bhattacharyya & Strohmayer 2005b). We argue that the Oct. 15 burst was ignited at a 
comparatively lower latitude, and hence v was greater than V for some time after equatorial 
belt ignition, while the Oct 19 burst and the bursts from 4U 1636-536 (that showed burst 
rise frequency evolution; Bhattacharyya & Strohmayer 2005b) were ignited at slightly higher 
latitudes. 

In their study of thermonuclear burning on rotating neutron stars Spitkovsky et al. 
(2002) derived the flame speed formula, z'flame = ( 9 h)^ 2 / ft n , where g , h, f and t n are 
the surface gravity, burning layer scale height, Coriolis parameter (/ = 20 spin cos 6), and 
nuclear burning timescale, respectively. An interesting question is whether the spin frequency 
dependence can be tested by comparing oscillation properties in bursts from different sources. 
Although one might expect such properties as burst rise times to reflect the flame speed, 
other time scales, such as the radiative diffusion time scale, are also relevant for rise times, 
and could mask the front speed effects. A perhaps more robust indication of the flame 
speed could be the time scale over which the oscillation amplitude drops during the early 
rising phase. Initial indications are that the amplitude drops more quickly in the Oct. 15 
burst from SAX J1808. 4-3658 than a burst from 4U 1636-536 that also shows frequency 
evolution during the rising phase (see Bhattacharyya & Strohmayer 2005b). At face value 
this is consistent with the spin frequency dependence from Spitkovsky et al. (2002), that 
is, other things being equal, flames propagate more slowly on more rapidly spinning stars. 
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However, this is still just suggestive, as it is difficult to be precise based on only several bursts. 
Nevertheless, it seems clear that further studies of burst rise frequency evolution from a larger 
sample of bursts may provide important statistical information about how the flame speed 
and other flame spreading properties can be influenced by stellar spin, compactness, and 
burst strength. Therefore, we emphasize that such studies provide a potentially powerful 
tool to understand thermonuclear flame spreading on neutron star surfaces under extreme 
physical conditions. 

This work was supported in part by NASA Guest Investigator grants. 
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Table 1. Frequency evolution model parameters a (with la error) for burst rise oscillations 

of three bursts from SAX J 1808. 4-3658. 


Burst date 


h 

fbi 

*>2 

v 2 

X 2 / dof 

Zf 

2002 Oct 15 

396 . 721 ^ 

29.64™ 

0.201JS 

-15.54itS 

23.67^ 1 8 f 

3.36/4 

152.16 

2002 Oct 18 

398.22l{j$ 

- 

0-3 lto.04 

19.2ll|39 

- 

8.44/6 

62.83 

2002 Oct 19 

399.64^014 

i-69±g;S 

- 


- 

6.46/7 

86.64 


a Frequency evolution model: u(t) = +v\t (for t < t^, and = z/(0)); v{t) = i/(f bl ) + z> 2 (t — 
tb x ) + - tbi) 2 (for t > tbi). 


b Fundamental power during burst rise. 
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Fig. 1. — Time evolution of different observed burst properties during the rise of the Oct 
15 burst from SAX J 1808.4-3658. Panel a gives the detected intensity (histogram), power 
contours (minimum and maximum power values are 16 and 51) from the dynamic power 
spectra (for 0.15 s duration at 0.01 s intervals), and the best fit model from Table 1. Panel b 
shows the rms amplitude of the oscillations. Here the horizontal lines give the binsize. Panel 
c gives the inferred radius (assuming 10 kpc source distance) of the source, while panel d 
gives the corresponding temperature. For panels 6, c, and d, persistent emission is subtracted 
and a deadtime correction is applied, and the error bars are la values. 
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Fig. 2. — Similar to panel a of Fig. 1 but for the Oct 18 (upper) and Oct 19 (lower) bursts. 
For the upper panel dynamic power spectra are calculated for 0.15 s duration at 0.01 s 
intervals, while for the lower panel these numbers are 0.2 s and 0.01 s. For the calculation of 
power contours, minimum and maximum values of power are 15 and 50 for the upper panel, 
and 17 and 111 for the lower panel. The dotted line in the upper panel gives the frequency 
evolution for which the total power (74) is higher than that (63) corresponding to the best 
fit model (solid line) from Table 1. 
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